Base Model

\[\begin{align} S' & = \mu N^\star - \lambda S - \mu S \\ E' & = \lambda S - \kappa E - \mu E \\ I_S' & = p \kappa E - (\gamma_S + \mu_{I_S} + \delta_H) I_S - \mu I_S \\ I_A' &= (1 - p) \kappa E - \gamma_A I_A - \mu I_A \\ H' &= \delta_H I_S - (\gamma_H + \mu_H) H - \mu H \\ R' & = \gamma_S I_S + \gamma_A I_A + \gamma_H H - \mu R \\ D' &= \mu_{I_S} I_S + \mu_H H \\ \lambda &:= \frac{\beta_A I_A + \beta_S I_S}{N^{\star}} \\ N^{\star}(t) &= S + E + I_S + I_A + H + R \end{align}\]

Stan ODE notation beginning model

Observation Model

\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t \rho \delta_E E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa: \sim \text{Gamma(10,50)} \end{align}\]

March 10 – March 30 parameters estimation

Postulated densities

knitr::opts_chunk$set(echo = FALSE)
library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

mu_ <- function(lambda, ylabel= "mu_h"){
    # We suppose mu_. as a exponential r.v.
    x <- seq(from = 0, to = 1, by=.001)
    mu_h_ <- dexp(x, rate = lambda)
    df <- data.frame(x=x, mu_h = mu_h_)
    mean_ <- 1/lambda
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
        p <- ggplot(data=df,
                aes(x = x,
                    y = mu_h)) + 
        xlab("x")+
        ylab(ylabel) +
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") + 
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5)+
        annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
              color="black")
    return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
    # We suppose mu_. as a exponential r.v.
    theta <- 1 / beta 
    x <- seq(from = 0, to = 1, by=.001)
    gamma_s_ <- dgamma(x, shape = alpha, 
                        scale= theta)
    df <- data.frame(x=x, gamma_s = gamma_s_)
    mean_ <- alpha * theta
    mean_time <- 1/mean_ 
    text <- paste(mean_time,"days")
    p <- ggplot(data=df,
                aes(x = x,
                    y = gamma_s)) + 
        geom_line(colour = "#000000") +
        geom_area(alpha = 0.6, fill = "lightgray") +
        geom_vline(xintercept = mean_, linetype="dotted", 
                color = "red", size=0.5) + 
        annotate(geom="text", x=1.30 * mean_, y=0, label = text,
              color="black")
        return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100,  ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50,  ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250,  ylabel= "gamma_H")

fig6 <- ggplotly(pi_mu_h)

## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6

fig7

fig8

fig9

fig10

Stan estimation

Lockdown-Vaccination-Treatment-Model

\[\begin{align} S'(t) &= \mu \bar{N} -\frac{\beta_S I_S+\beta_AI_A}{\bar{N}}S -(\mu+\lambda_V)S +\delta_V V \\ E'(t) &= \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} S + \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}}V - (\mu+\delta_E) E \nonumber \\ I'_S(t) & = p \delta_E E -(\mu + \mu_S + \alpha_S + \lambda_T) I_S +(1 - q) \alpha_T T \\ I'_A(t) & = (1-p) \delta_E E - (\mu + \mu_A + \alpha_A) I_A \\ R'(t) & = \alpha_S I_S + \alpha_A I_A + q \alpha_T T - \mu R \\ D'(t)& = \mu_S I_S + \mu_A I_A \\ V'(t)& = \lambda_V S - \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}} V - (\mu + \delta_V) V \\ T'(t) &= \lambda_T I_S -(\mu + \alpha_T) T \\ \bar{N}(t) &= S(t) + E(t) + I_S(t) + I_A(t) + R(t) + V(t) + T(t) \end{align}\]

Model from March 24 July 14

March 24 –Jun 01 parameters estimation

Stan ODE notation

Stan MODEL

## functions {
##     real[] seirvt(real t,  // time
##              real[] y,
##              real[] theta,
##              real[] x_r,
##              int[] x_i) {
## //
##     real dy_dt[8];
##     real s = y[1];
##     real e = y[2];
##     real i_s = y[3];
##     real i_a = y[4];
##     real h = y[5];
##     real r = y[6];
##     real d = y[7];
##     real n_star;
## //
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real kappa = theta[3];
##     real p = theta[4];
##     real delta_h = theta[5];
##     real mu_i_s = theta[6];
##     real mu_h = theta[7];
##     real gamma_s = theta[8];
##     real gamma_a = theta[9];
##     real gamma_h = theta[10];
## //
##     real force_infection;
##     real mu = 3.913894e-05;
##     n_star = s + e + i_s + i_a + h + r;
## 
##     force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
##     dy_dt[1] = mu * n_star - force_infection * s - mu * s;
##     dy_dt[2] = force_infection * s - (mu + kappa) * e;
##     dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
##     dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
##     dy_dt[5] = delta_h * i_s  - (gamma_h + mu_h + mu) * h;
##     dy_dt[6] = gamma_s * i_s + gamma_a * i_a  + gamma_h * h - mu * r;
##     dy_dt[7] = mu_i_s * i_s + mu_h * h;
##     dy_dt[8] = p * kappa * i_s;
##     return dy_dt;
##   }
## }
## data {
##     int<lower = 1> n_obs;       // number of days observed
##     int<lower = 1> n_theta;     // number of model parameters
##     int<lower = 1> n_difeq;     // number of differential equations
##     int<lower = 1> n_pop;       // population
##     int y[n_obs];               // data, total number of infected
##     real t0;                    // initial time point (zero)
##     real ts[n_obs];             // time points observed
## }
## 
## transformed data {
##   real x_r[0];
##   int x_i[0];
## }
## parameters {
##     real <lower = 5e-6> theta[n_theta];
##     // real <lower = 0, upper = 1> E0;
##     // real <lower = 0, upper=1> IA0;
##     real <lower = 0, upper = n_pop> E0;
##     real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
##     real y_hat[n_obs, n_difeq]; // solution from the ODE solver
##     real y_init[n_difeq];
##     // initial conditions
##     // real IS0 = 8.297217e-06;
##     real IS0 = 22.0;
##     real H0 = 0.0;
##     real R0 = 0.0;
##     real D0 = 0.0;
##     real CIS0 = 74.0;
##     y_init[1] = n_pop - (IS0 + IA0 + E0);
##     y_init[2] = E0;
##     y_init[3] = IS0;
##     y_init[4] = IA0;
##     y_init[5] = H0;
##     y_init[6] = R0;
##     y_init[7] = D0;
##     y_init[8] = CIS0;
##     y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
##     real lambda[n_obs];             // poisson parameter
##                                     // priors
##     theta[1] ~ normal(1.0, 0.1);    // beta_s
##     theta[2] ~ normal(1.0, 0.1);    // beta_a
##     theta[3] ~ gamma(10, 40);       // kappa
##     theta[4] ~ uniform(0.1, 0.5);   // p
##     theta[5] ~ gamma(10, 40);       // delta_h
##     theta[6] ~ exponential(33.33333);   // mu_i_s
##     theta[7] ~ exponential(25);   // mu_h
##     theta[8] ~ gamma(10, 100);      // gamma_s
##     theta[9] ~ gamma(10, 50);       // gamma_a
##     theta[10] ~ gamma(10, 250);      // gamma_h
## //
##     // E0 ~ uniform(1.121246e-07, 3.363737e-05);
##     // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
##     E0 ~ uniform(74, 2100);
##     IA0 ~ uniform(74, 2100);
## //    likelihood
##     for (i in 1:n_obs){
##        //lambda[i] = y_hat[i, 8] * n_pop;
##         lambda[i] = y_hat[i, 8];
##     }
##     y ~ poisson(lambda);
## }
## generated quantities {
##     real R_0;      // Basic reproduction number
##     real mu = 3.913894e-05;
##     real beta_s = theta[1];
##     real beta_a = theta[2];
##     real p = theta[4];
##     real gamma_s = theta[4];
##     real gamma_a = theta[5];
##     real kappa = theta[6];
##     R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
##         + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }

Observation Model

\[\begin{align} Y_t & \sim \text{Poisson}(\lambda_t), \qquad \\ \lambda_t & = \int_0^t p \kappa E \\ & p \sim \text{Uniform(0.3 0.8)} \\ & \kappa \sim \text{Gamma(10, 52)} \end{align}\]